% By Martin M. Andreasen, April 22 2010
% This function simulates the model when solved up to third order. The purning scheme is NOT used.
% These codes are faster for a third order approximation because we use the kronecker product to eleminate loops
function [Y_sim,X_sim] = simulate_3rd_kron_NoPruning(gx,gxx,gxxx,hx,hxx,hxxx,gss,gssx,gsss,hss,hssx,hsss,sig,eta,shocks,order_app)

% Some indices 
ny      = size(gx,1);        %Number of output variables
nx      = size(hx,1);        %Number of state variables
ne      = size(eta,2);       %Number of structural shocks
num_sim = size(shocks,2);    %Number of simulations

% Allocating memory
Y_sim = zeros(ny,num_sim);
X_sim = zeros(nx,num_sim);

% The simulation
xt    = zeros(nx,1);    
AA    = kron(xt,xt);

% Defining matrices
Hxx  = 1/2*reshape(hxx,nx,nx^2);
Hxxx = 1/6*reshape(hxxx,nx,nx^3);
Gxx  = 1/2*reshape(gxx,ny,nx^2);   
Gxxx = 1/6*reshape(gxxx,ny,nx^3);   

if order_app == 1
   for t=1:num_sim
      % The state variables
      xt_p  = hx*xt + sig*eta*shocks(:,t);
      X_sim(:,t) = xt_p;

      % The controls
      Y_sim(:,t) = gx*xt_p;
    
      % Updating xt
      xt  = xt_p;
   end
elseif order_app == 2
   for t=1:num_sim
      % The state variables
      xt_p  = hx*xt + sig*eta*shocks(:,t)...
            + Hxx*AA + 1/2*sig^2*hss;
      X_sim(:,t) = xt_p;

      % The controls
      AA = kron(xt_p,xt_p);
      Y_sim(:,t) = gx*xt_p + Gxx*AA + 1/2*sig^2*gss;
    
      % Updating xt
      xt  = xt_p;
   end
elseif order_app == 3
   for t=1:num_sim
      % The state variables
      xt_p  = hx*xt + sig*eta*shocks(:,t)...
            + Hxx*AA + 1/2*sig^2*hss...
            + Hxxx*kron(xt,AA)+ 3/6*hssx*sig^2*xt + 1/6*sig^3*hsss;
      X_sim(:,t) = xt_p;    
      
      % The controls
      AA    = kron(xt_p,xt_p);
      Y_sim(:,t) = gx*xt_p + Gxx*AA + 1/2*sig^2*gss + ...
                   Gxxx*kron(xt_p,AA) + 3/6*gssx*sig^2*xt_p + 1/6*sig^3*gsss;
        
      % Updating xt
      xt  = xt_p;
   end
end
